Polynomial splines interpolating prime series 

L. Alexandrov* D. B. Baranov* and P. T. Yotov^ 

Abstract: Differentiable real function reproducing primes up to a given num- 
ber and having a differentiable inverse function is constructed. This inverse 
O ' function is compared with the Riemann-Von Mangoldt exact expression for the 

^ ■ number of primes not exceeding a given value. Software for computation of 

^ . the direct and inverse functions and their derivatives is developed. Examples 

Q I of approximate solution of Diophantine equations on the primes are given. 

^ ■ 1. Introduction 

(N 

This article introduces real functions reproducing the values of mutually inverse arith- 
f-H I metic functions p{n) : N ^ P {prime p{n) at a number n) and p~^{p) : P ^ N {number 

^ ■ n of the prime p{n)). 

"5 ! The found functions are employed to create subroutines for computation of '^^^'^ 

on 1 < a; < cxd, and p'^^{x)^ dx on 2 < x < oo. 

' The above-noted programs can be used for a numerical solution of different problems 



> 



(N 
(N 



X 



on the set of primes P, including approximate solution of Diophantine equations on P. 



Tj- . The idea consists in establishing a differentiable function which would include the 



values of primes and which would allow one to construct an inverse function p [x) by a 
Newton method. More precisely, the sought function 1 < x < oo should satisfy the 
^ I following conditions: 

■ a) p{x) reproduces primes; 

\ b) there exists a positive derivative 

^ ! c) there exists an inverse function p {x) : [2, oo) [1, oo). 

j>I I As is known, there are no one-variable polynomials which can produce all the primes, 

or primes only. However, this article shows that there exist polynomial splines reproducing 
primes in series (along with the continuation of the prime series) and satisfying in addition 
the conditions b) and c). 

A spline formed by polynomials with integer coefficients will be called the arithmetic 
spline. 

This article discusses two candidates for the arithmetic splines, cubic one and parabolic. 
These splines do not approximate a prime series. Primes are implanted in the structure 
of the splines, which ensures that they are exactly reproducible. Such splines lead to ex- 
plicit soluble systems of linear equations whose coefficients represent arithmetic functions 
themselves. 

The inverse function p~^{x) constitutes a differentiable analogue of the number- 
theoretic function 7r(x) = 'Yl,^ which is comparable with the Riemann exact expression 

for 7r(x) through the zeros of the C^function( pp, page 34). 
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2. Cubic spline 

Consider the sphne 



( x + 1, l<x < 1.5, 

Ci{x), i — 0.5 < X < i + 0.5, i = 2,3,..., 
Ci{x) = Ci+i{x) , x = i + 0.5, i = l,2,..., 



with 



Ci{x) = 2 {ai{x - i - 0.5)2 ^ ^.(^ _ . _ ^ + - «)- 

-2p(i)(a;-i-0.5). 
Exact reproducibihty of the primes follows from the identity 

Q(i)=p(i), i = l,2,.... (1) 

At the points of sewing the spline should also obey the identity 

Q(i + 0.5) = ^(p(i)+p(i + l)), i = l,2..., (2) 

which brings into the spline additional information of the prime series behaviour. 
There exists an unique cubic spline of the kind Scubi^) with the coefficients 



ai^^ip{i + l)-p{i-l))-l, 



bi ^p{i + 1) -p{i) - 1, 



i = 2,3... 



This spline can be considered only as almost-arithmetic. The coefficients 7^ and 5i in 
Ci{x) — aiX^ + j3ix'^ + ^iX + 5i appear for some i in the form q + 1/2, q & N. 
The negative value of the discriminant 

di = 4{p{t)f - 4{p{t - 1) + p{t + l))p{i) + ^{pit - 1))2 + ^{p{t + 1))2+ 

+ ^p{i - l)p{i + 1) + 3 



gives positivity of the derivative 

rln-( 1" 1 

= 2{2ai{x - i - 0.5) + bi){x - i) + 2ai{x - i - 0.5f + 2bi{x - i - 0.5)- 



dx 

-p{i)+p{i + l). 

The inequality di < leads to the following condition for prime triplets 
p{i - 1), p{i), p(i + l): 

_ 1) + 1)) - i < p(^) < - 1) +p(t + 1)) + 1 Vt;, (3) 




(4) 



where 

U = SMt + 1) - Pit - \)f - 4) > 0, ^ = 1, 2, . 
Condition Q can be violated in the cases 

- 1) = Ai, p(z + 1) > A2; 

p(z) - 1) > A2, p(i + 1) = Ai, 

with Ai = 2, A2 = 28. 
Among the first 1000 primes only 5 triplets violate the rule Q: 

(2969, 2971, 2999), 
(2971, 2999, 3001), 

(3271, 3299, 3301), (5) 
(6917, 6947, 6949), 
(7757, 7759, 7789). 

Except these triplets (including twin pairs), condition is violated by triplets of the 
kind (0} at the following values for Ai and A2 : 

Ai = 4, A2 = 56, 
Ai = 6, A2 = 84, 
Ai = 8, A2 = 114, 
and all that. 
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Figure 2: Derivatives S'^„^(a;)(quadrics) and ^^^^^(x) (segments) 

Despite the cases where condition (jH} is violated, sphne Scub{x) is convenient for cre- 
ating subroutines p{x), '^^^j^'^ and p^^{x), since the inverse function p~^[x) exists in the 
neighborhood of each prime number. 

3. Parabolic spline 

Given the following pairs of parabolas 



"5 " '-'5 

i = 2,3, 



' qlix)^ i — 0.5<x<i, i = 2,3, 
i < X < i + 0.5, 

^[■^j - Vii^; I 

X = 2, 3, ... ; internal sewing. 



ql[x) = q-[x) 



dQiix) ^ dql{x) 

dec doc 



with 



qi{x] 

ql{x 
ai = 



2aj_i(a; — iY + x — i + p{i) 
ai{x — i — 0.5) 
p{i + 1) -p{i) - 1. 



2ai{x - % - 0.5)2 ^ (2ai + l)(x - % - 0.5) + + 
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The parabolic spline 



' X + 1, 1 < X < 1.5, initial polynomial, 
qi{x), i — 0.5 < X < i + 0.5, i = 2,3, . . . 

g.(x)=g.+i(x) I a: = ^ + 0.5, i- 
dqijx) _ dqi+i{x) [ ' external sewing 



1,2,...; 



dx 



dx 



solves the problem better than the spline Scub- It has the following properties: 
1) identities analogous to ((1} and (jSj) are applicable 

qi{i) =p{i). 



(6) 



g,(^ + 0.5) = -(p(2)+p(z + l)), z = l,2. 



2) the derivatives 



dq\{x) 

dx 

dql{x) 

dx 



4aj_i(i - x) + 1, 



AaAx — i) + 1 



(7) 



take a minimal value +1 at the points of internal sewing (they are the points of inter- 
polation to the spline) and maximal values at the points of external sewing, where they 
coincide with the derivative of the spline Scuh- 
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The positive values of the derivatives show that the sphne Squad{x) monotonically 
increases on the semi- axis [1, oo). 

There exists a function S~^^^{x), inverse to the function Squad{x), determined on the 
axis [2, oo) and thus the sphne Sguad fulfiUs the conditions a), b), c). Moreover, the 
function S~^^^{x) is differentiable on (2, oo). 

The sphne Squad is arithmetic because prime number polynomials 

ql{x) = a\x'^ + (3[x + 7' and ql{x) = + (3lx + 7[ 
hold integer coefficients (see Table 1): 

a\ = -2aj„i; = 4iaj_i + 1; 'yl = -2i^aj_i + p{i) - i; 

a[ = 2aj; /3[ = -Aiai + 1; 7[ = 2i^ai + p{i) - i. 

The joint satisfiability of identities ((Tj) reinforces the hypothesis that the spline Squad 
is an unique arithmetic spline satisfying the conditions a), b) and c). 



Table 1 {dl = iPir - AcH. dl = (AO' - 4al7r) 



i 


p{i) 




PI 
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<: -y^ 


dl 




2 


3 





1 


1 


1 


2 


-7 


9 


-23 


3 


5 


-2 


13 


-16 


41 


2 


-11 


20 


-39 


4 


7 


-2 


17 


-29 


57 


6 


-47 


99 


-167 


5 


11 


-6 


61 


-144 


265 


2 


-19 


56 


-87 


6 


13 


-2 


25 


-65 


105 


6 


-71 


223 


-311 


7 


17 


-6 


85 


-284 


409 


2 


-27 


108 


-135 


8 


19 


-2 


33 


-117 


153 


6 


-95 


395 


-455 


9 


23 


-6 


109 


-472 


553 


10 


-179 


824 


-919 


10 


29 


-10 


201 


-981 


1161 


2 


-39 


219 


-231 


11 


31 


-2 


45 


-222 


249 


10 


-219 


1230 


-1239 


12 


37 


-10 


241 


-1415 


1481 


6 


-143 


889 


-887 


13 


41 


-6 


157 


-986 


985 


2 


-51 


366 


-327 


14 


43 


-2 


57 


-363 


345 


6 


-167 


1205 


-1031 


15 


47 


-6 


181 


-1318 


1129 


10 


-299 


2282 


-1879 


16 


53 


-10 


321 


-2523 


2121 


10 


-319 


2597 


-2119 


17 


59 


-10 


341 


-2848 


2361 


2 


-67 


620 


-471 


18 


61 


-2 


73 


-605 


489 


10 


-359 


3283 


-2439 


19 


67 


-10 


381 


-3562 


2681 


6 


-227 


2214 


-1607 


20 


71 


-6 


241 


-2349 


1705 


2 


-79 


851 


-567 


21 


73 


-2 


85 


-830 


585 


10 


-419 


4462 


-2919 



It should be noted here that the coefficient in ql{x) and g[(x) represents a basic 
arithmetic function - the number of composite numbers in the interval {p{i),p{i + 1)). 

Figures 1-3 present a comparison of the splines Scub and Squad- The interval's [428, 432] 
image (argument to the functions p{x) and dp{x)) contains the first pair of triplets (0) 
violative the positivity of the derivative 5*^^,6 (x). 
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Figure 3 shows intervals where no inverse function p~^{x) for the sphne Scub exists. 
Figure 2 illustrates the properties of the derivatives S'^^f^ and 3'^^^^^: reaching the min- 
imal and maximal values and the equalities of derivatives at the sewing points. 

4. Inverse parabolic spline S~^^^{x) and its derivative 

The pairs of functions 



where 



ti{x) 



tKx), MO < X < ^«±#±1) , i = 2,3, 



dx dx 



a; =p(2),p(3),..., 



1 - \ h\ 



t'i{x)=i+ 4^.:^ ,(4 = x-l), t\{x) = i + 
h\ = 8ai_i(p(i) - x) + 1, %ai{x - p{i)) + 1, 

ai =p{i + 1) -p{i) - 1, 



6r-i 



determine the inverse spline 



( x-1, 2<x < 2.5, 
p{i - 1) +p{i) 



(iti(3;) ^ dtj+i {x) 
dx dx 



<x< 



p{l)+p{i + l) 



, i — 2,3, ... , 



p{i) +p{i + 1) 



,z = l,2,.... 



The derivative of the inverse spline is as follows 



dSJadi^) 



quad\ 

dx 



= < 



1, 2<x<2.5, 

,a-i/2 p{i - 1) +p{i) 



2 < a; < i = 2, 3, . . . , 

, p(i)<.<£ei±|(i+i),,=2,3,.... 

5. About subroutines p(x), and ^^^J^-^ . 



Both splines ^c^b and ^'^^ad were employed to create Fortran functions p{x), 
dp{x) :— P-{x)i P-newt{x) {p-{x) and p^newt{x) denote p~^{x)) and 

, / s dp'^ix) 
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59450 - 



59400 



59350 - 



59300 



59250 



59200 - 



59150 



5980 5985 5990 5995 6000 6005 



6010 



Figure 4: Sewing the spline Sguad{x) with asymptote p{x) at x=6000. 

For convenience in apphcations the functions p{x) and dp{x) have been extended to 
+00 by the asymptote (j2], page 140) 

, Inlnx — 2 (lnlnx)^/2 — 31nlnx + 5.5 
p[x) = X In X + In In a; H : — 1 | . 



Inx 



(Inx) 



To obtain primes, two alternative programs have been employed - subroutine eratos- 
thenes(n) and subroutine primes(n). The former accomplishes this by generating primes 
up to a given value n, whereas the latter achieves the result by reading n primes from 
given 6 column file named primes. In these programs, the spline and its derivative are 
automatically sewed with asymptote (jH)) and its derivative. 

For purposes of building the inverse function p~^{x) in the programs p^newt{x) and 
P-{x), two different approaches have been used. 

The program p_newt{x) is based on an autoregularized variant of the Newton method 
(|S],page 43). 



I/O, £0 > 0, Vk+i = Vk - 



X 



dpiVk) + Sk ' 
£fc = 2 [\/{dp{yk)f + ^N\p{yk) - x\ - dp{yk 



k = 0,1,2, 



(9) 



N = {el + eodp{yo))/\piyo) - x\, 
where several combinations of the initial value yo = li{x) and initial regularizator Eq (these 
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combinations can be seen at the beginning of the p-newt{x) program's body) ensure a 
construction of the inverse function p~^{x) on the interval [2, 10*]. 
The way of setting the initial approximation yo within the values 

2 

= V ^/(x^/") (see iU, page 35), (10) 
^-^ n 

n=l 

where f{x) = li{x), fi{n) is if n is divisible by a prime square, 1 if n is a product 
of an even number of distinct primes, and —1 if n is a product of an odd number of 
distinct primes was checked: the conclusion is that the method Q and = li{x) are an 
acceptable combination. 

Figure 4 shows the sewing of the function p{x) with the asymptote (jH)). 

The programs p_(x) and dp-{x) are based on the application of the inverse spline 
^quad derivative. In these programs, an algorithm has been applied by virtue of 

which the needed pairs ti{x) and bi{x) are established by approximating ix = [li{x)\ (see 
the beginning of the programmes p-{x) and dp-{x)). It is worth mentioning that this 
algorithm works in the case where tt{x) < li{x), as well as in the case where 7r(a;) > li{x), 
i.e., the algorithm does not depend on the knowledge the minimal value of x for which 
the difference li{x) — 7i{x) changes the sign. 

The programs p{x), dp{x), p-newt{x), p^{x) and dp^{x) have been realized in the 
Fortran90. These programs form the package named pp„/90[l, oo), available in Appendix 
1. They employ the natural extension of the function to — oo, based on the values of the 
initial polynomial x + 1. 

The package pp_/90[l, oo) is immediately applicable to Compaque- and MS- Fortran 
and is facile transportable to other Fortran versions. 

Fortran functions p{x) and p~^{x) are as easily applicable as the intrinsic functions 
sin{x) and exp{x). 

The end of Appendix 1 contains a program called test^pp^ which has been used to 
compute all the tables supporting the graphics in this article and which serves as an il- 
lustration for application of the functions p{x) and p~^{x). 

6. Possible applications of the functions p{x) and p~^{x). 

6.1. Functions p{x) and p~^{x) can be used to introduce new functions sinp(x), 
cosp{x), tanp(x), e'^^^^ and Inp(x), apphcable when the specific character of nonasymp- 
totic prime number distribution is necessary to be accounted for (see, e.g., H]). In par- 
ticular, the functions sinp(x) and cosp{x) can be used for creation of a prime number 
harmonic analysis. 

6.2. Diophantine equations can be solved within the following approximate method: 
to the Diophantine equation 

/i(xi,...,xO = (11) 
one adds a new equation, from reals-to-integers equation 

/^(xi, ...,Xn) := sin^(7rxi) + sin^(7rx2) + ■ ■ ■ + sin^(7rx„) = 0, (12) 
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or one adds from reals-to-primes equation 

f2ixi,...,Xn) := sm2(7rp"^(a;i)) + sin^ (np''^ {X2)) H h sin^(7rp"^(a;„)) = 0. (13) 

By solving either system (fTT|). (fT^ or (fTTj) . (fT^ one can find solutions to (fTT|) as real 
approximations to the natural or the prime numbers. 

The above systems can be solved by methods working in the case of degeneration of 
the derivative at the solution (see, e.g., jS]). 

Here two examples of solving such systems by means of autoregularized iterative pro- 
cesses {rgn) [H] are presented; rgfn-pro cesses are combined with both the svd-method [7j 
and the adaptive scaling |H]. 

In this case the program afxy [1^ is used to find all solutions of the systems (fTT| . (fT2|l 
and (jllj) . (fT!?|l in a given definition domain. 

Let the problem (|TT|l . (|T!^ be written as 

Fx = 0, (14) 

where 

Fx = f'^{x){fx-y), fx = [A (a;), ^(x)]^, f : D f G ^ , x e ,y e R^, 

and is an open convex domain in R^. 

The linear problem at the fcth iteration of r(7n-process is of the kind 



(/"^(a;')/'(x'^) + £fc/)(x'=+i - a;^) = -Fx\ (15) 



where 



1 



n = ||/'^(x'=)/'(x'^)|U, Pk = \\Fx^\U, iV = (£0 + £oro)/po. 

For simplicity in equality (jl5p and in the expressions for t^, pk and £^1, the scaling oper- 
ators are neglected. Just for the problem (jl5p the sf d-method is in use. 

For the purpose of finding all solutions of equation ()14|) in Df, the program afxy 
realizes an algorithm in which the vector Fxk is multiple factored by local extractors of 
the kind 

eJx.x^'^) = (1 



where x^'^^ E Df is the rth solution of equation (fT^ . The transformed problem 

Fr^x := lf[er{x,x^''^)\ Fx = 0, r* > 1 (16) 



vr=l 



is multiple solvable by the program afxy. 

For each new problem (fT^ afxy realizes different rgn iterative processes: with different 
guesses Xq (randomly formed by the initially given xq and by peculiarities of the domain 
Df) and with different initial regularizators Eq (from the preset table of regularizators) . 

Example 1. Solution of the equation xl + X2 = x"^ + 1 over primes. 
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Consider problem ()14|) with 

( fi{x) ■.= xl + xl-xl-l = 0, 

[ f2{x) := sin^{np ^{xij) + sin^(7rp ^(0:2)) H h sin^(7rp = 0, 

n = 3,m = 2,Df = [2, 100] x [2, 100] x [2, 100], x = (xi, Xa, x^f , y = [0, 0]^. 

Run-time section of Application 2 contains a subroutine fxy{m, n, np, neq, f, x, pp, df, yr), 

dfix) 

where the equation fx = and derivatives (i = 1, 2, 3, j = 1, 2) are coded. Here, in 

the pre-exe section of the program afxy, the main controls m,n,np, f,lsmh,mqh,nsolh, 
x^i\xf^ and Df are given as well. 

If in the pre-exe section some main control is not prescribed, then afxy switches to an 
interactive mode and demands from display an adjustment to the value of this control. 

The solution (xi,X2,X3) of the equation /i(x) = ([ni, page 35) is thought to be a 
quasi- Pythagorean prime triplet (the equation x1 -\- x^ = x| has no solutions on primes). 
There exist two series of natural numbers satisfying equation fi{x) = : 

(2n + 1, ra^ + n - 1, ra^ + n + 1) , (TUl, page 35, 

{2n{An + 1), IQn^ - 1, IQn^ + 2n) , [lOj, page 36. 

The present state of the program afxy can only produce up to 20 (i.e. r* < 20) 
extractions and can find only 20 quasi- Pythagorean prime triplets in Df{see FOUND SO- 
LUTIONS in Application 2). 

Example 2. Solution of the equation x^ + X2 = x| + 1 over twins. 

Consider problem (fT^ with 

/i(x) := X? + x| -x| - 1 = 0, 

fx = < /2(x) := sin^(7rp"^(xi)) + sin^(7rp"^(x2)) H h sin^(7rp~^(x„)) = 0, 

, fsix) := X3 - xi - 2 = 0, 

n = 3,m = 3,Df = [2, 100] x [2, 100] x [2, 100], and y = [0, 0, 0]^. 

The needed subroutine fxy for this example is presented in Application 3. 

Application 3 shows that the program afxy finds all 5 quasi- Pythagorean prime triplets 
containing twin pairs in the domain Df. 

6.3. In Figure 5, the inverse function p~^{x) is compared with the step-function 7r(x) 
and with the Riemann-Von Mangoldt continuous step-function 7ri^(x) which results from 
(fTUj) by means of the substitution 



00 

Uix)-f^U{ePn^^^)+ I , -ln2, 

J [a —lama 
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where {pn} are the complex zeros of the equation 



in the form p„ = ^ ± itji, n 




(17) 



In our case, counting function 7i{x) is expressed by the inverse function p ^{x) as the 
formula 



Comparing the two ways, the function p^^{x) proves more convenient for application 
than the Riemann-Von Mangoldt function 7r/j(x); moreover, p{x) and p~^{x) are appli- 
cable just now, not waiting for the final representation of the function vr/j(x) through the 
zeros Pn- 

6. 4. For purposes of investigation the nonasymptotic behaviour of primes, instead of 
the functions p{x) and p^^{x), one can resort the functions: 

A{x) = p{x) - p{x) - (p(xo) - pixo)), X e [xo, xo + e{xo)) 

-a local variance of the function p{x), where e{xo) is relatively small in comparison with 
xq, and 



-a local variance of the function p^^{x). 

Here p{x) is the asymptote (jH)), and R{x) is the Riemann simplified formula ()1U|) for tt{x). 
Figures 7-9 serve as examples to the behaviuor of A{x) and B{x). 



7l{x) = \j) (x)J. 



(18) 



B{x) = p ^{x) - R{x) - {p ^{xo) - -R(xo)), X E [xo, xq + e{xo)) 
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Figure 6: Function p{x). The interval [154.78, 168.2] is chosen so as the definition interval 
[900, 1000] of the function p~^{x) coincides with the definition interval of the function 
number variance of the zeros pn from Figure 4 [T^, p. 406. 





A 

A : 

158 J 62.- 




\ 156 

V 




164 166 168 



Figure 7: Function A{x); note a qualitative intimacy with the function number variance 
of the zeros pn from Figure 4 [TT], p. 406, as well as the same number of peaks equal 
to 14 in the given interval; the function A{x) covers 14 successive primes: 907, 911, 919, 
929, 937, 941, 947,953, 967, 971, 977, 983, 991 and 997. 
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Figure 8: Function B{x). Inverse function to the function A{x) from Figure 7. 




Figure 9: Function B{x) on the enlarged interval [900,1150]. How many peaks does the 
function number variance of the zeros pn have on the interval [1000, 1150]? Do we really 
have no 17 peaks? 
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Application 1. Programm package ]9]9_/90[l, oo) 



function p(x) ! 1 <= x < infinity 
implicit real*8(a-h,o-z) 
common/ wn/wn 

if(x <= 2.d0) then 
p = x+l.dO; return 

endif 

if(x > 2.d0 .and. x <= wn) then 

ix=floor(x+0.5d0) ; xx=df loat (ix) ; p= sq(ix,x,xx) ; return 
endif 

if (x > wn) then 

p= r (x) ; return 
endif 
contains 

function sq(ix,x,xx) ! left-right quadric pair 

common/protiarithmi/qdOOOOOOO) 

if (x <= xx) then 

sq=-2 . dO* (q(ix) -q(ix-l) -1 . dO) * (x-xx) **2+x-xx+q(ix) 
else 

sq=2 . dO* (q(ix+l) -q(ix) -1 . dO) * (x-xx-0 . 5d0) **2+ & 
(2 . dO* (q(ix+l) -q(ix) ) -1 . dO) * (x-xx-0 . 5d0) + (q(ix) +q(ix+l) ) /2 . dO 
endif 

end function sq 
end function p 

function dp(x) ! 1 <= x < infinity ; derivative of p(x) 

implicit real*8(a-h,o-z) 

common/wn/wn 

if(x <= 2.d0) then 
dp = 1 . dO ; return 

endif 

if(x > 2.d0 .and. x <= wn) then 

ix=f loor (x+0 . 5d0) ; xx=df loat (ix) ; dp= dsq(ix,x,xx) ; return 
endif 

if (x > wn) then 

dp=dr(x) ; return 
endif 
contains 

function dsq(ix,x,xx) ! derivative of left-right quadrics 
common/protiarithmi/q( 10000000) 
if (x <= xx) then 

dsq=-4 . dO* (q(ix) -q(ix-l) -1 . dO) * (x-xx) +1 . dO 
else 

dsq=4 . dO* (q(ix+l) -q(ix) -1 . dO) * (x-xx-0 . 5d0) + & 
2 . dO* (q(ix+l) -q(ix) ) -1 . dO 

endif 

end function dsq 
end function dp 
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function p_newt(u) ! 2 <= x < infinity 

implicit real*8(a-h,o-z) 

parameter (xyO=l . d3 , xepsl=58 . d3 , epsl=3 . dO , xeps2=17 . d5 , 
eps2=300 , eps3=700 . dO , ytol=l . d-11 ,dytol=l . d-11 ,ktol=1000) 
common/rcorr/rc , drc/kk/k/ rw/ r 
drcw=drc ; drc=0 . dO 
x=u ; xw=x 
if(x >= 3.d0) then 

if (x < xyO) then 
y=x/dlog(x) 

else 

y=dlogintegral (x) 

endif 
else 

p_newt=x-l .dO 

return 
endif 

if( X < xepsl) epsO=epsl 
if (xepsl <= X < xeps2) eps0=eps2 
if (xeps2 <= x) eps0=eps3 
dtO=dp (y) ; rO=dabs (p (y) -x) 
en= (eps0**2+eps0*dabs (dtO) ) /rO 
r=0.dO; k=0; small=l.d200 
k=k+l 

yy=y; rr=r; t=p(yy)-x; dt=dp(yy) ; r=dabs(t) Icurrant value of 
eps=0 . 5d0* (dsqrt (dt**2+4 . dO*en*r) -dabs (dt) ) ! autoregularizator 
y=yy-t/ (dt+eps) ! autoregularized Newtonian iterator 

dif =dabs (y-yy) 
dr=dabs (r-rr) 
if(r >= ytol .and. k <= ktol .and. dr > dytol) then 

if (dif <= small) then 

small=dif ; ybest=y; rr=r 

endif 

goto 1 
else 

p_newt= ybest 
endif 
drc=drcw 
return 
end 



initial value of 
initial value of 
autoregularizator 
constant for 

autoregularizator formula 
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function p_(x) ! 2 <= x < infinity 
implicit real*8(a-h,o-z) 

common/protiarithmi/qC 10000000) /wn/wn/ck/vkoch/nn/nn 
1 if(x >= 3.d0 .and. x <= q(nn)) then 

if(x <= l.d3) ixx=f loor(x/dlog(x) ) 

if(x > l.dS) ixx=f loor(dlogintegral(x)-vkoch*dsqrt(x)*dlog(x)) 
if(q(ixx) == O.dO) qq=r (df loat(ixx) ) 
if(q(ixx) /= O.dO) qq=q(ixx) 
if (qq <= x) then 

is=-l; si=-l.dO 
else 

is= 1; si= l.dO 
endif 
else 

if(x < 3.d0) then 

p_=x-l.dO; return 
endif 
endif 

do i=l, nn 

ii=ixx-is*i 

if (si*(q(ii)-x) <= O.dO) then 
ix=ii; goto 3; endif 
enddo 

3 if (dabs(q(ix)-x) >= dabs (q(ix+is) -x) ) ix=ix+is 

xx=q(ix) ; xi=df loat (ix) ; p_= y (ix,xi ,x,xx) 
if (vkoch /= O.dO) goto 1 
contains 

function y(ix,xi ,x,xx) ! left-right inverse quadric pair 
implicit real*8(a-h, o-z) 
common/protiarithmi/qdOOOOOOO) /ck/vkoch 

vkoch=0 . dO 
if(x <= xx) then 

a=q(ix)-q(ix-l)-l.dO; b=8.d0*a*(q(ix)-x)+l .dO 
if(b <= O.dO) then 

vkoch=l .dO 
else 

if (a == O.dO) then 

y=x-l.dO; else; y=xi+(l .dO-dsqrt (b) ) / (4 . dO*a) ; endif 

endif 
else 

a=q(ix+l)-q(ix)-l.dO; b=8 .dO*a* (x-q(ix) )+l .dO 
if(b < O.dO) then 

vkoch=l .dO 
else 

y=xi+ (dsqrt (b) -1 . dO) / (4 . dO*a) 
endif 
endif 
end function y 
end function p_ 
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function dp_(x) ! 2 <= x < infinity; inverse of dp(x) 
implicit real*8(a-h,o-z) 

common/protiarithmi/qC 10000000) /wn/wn/ck/vkoch/nn/nn 
if(x >= 3.d0 .and. x <= q(nn)) then 
if(x <= l.d3) ixx=f loor(x/dlog(x) ) 

if(x > l.dS) ixx=f loor(dlogintegral(x)-vkoch*dsqrt(x)*dlog(x)) 
if(q(ixx) == O.dO) qq=r (df loat(ixx) ) 
if(q(ixx) /= O.dO) qq=q(ixx) 
if (qq <= x) then 

is=-l; si=-l.dO 
else 

is= 1; si= l.dO 
endif 
else 

if(x < 3.d0) then 

dp_=l.dO 

return 
endif 
endif 

do i=l, nn 

ii=ixx-is*i 

if (si*(q(ii)-x) <= O.dO) then 
ix=ii; goto 3; endif 
enddo 

if (dabs(q(ix)-x) >= dabs (q(ix+is) -x) ) ix=ix+is 
xx=q(ix) ; xi=df loat (ix) ; dp_= dy(ix,x,xx) 
if (vkoch /= O.dO) goto 1 
contains 

function dy(ix,x,xx) ! left-right inverse quadric pair 

coinmon/protiarithmi/q(10000000)/ck/vkoch 

vkoch=0 . dO 

if(x <= xx) then 

a=q(ix) -q(ix-l) -1 . dO ; b=l . dO+8 . d0*a* (q(ix) -x) 

if(b <= O.dO) then 
vkoch=l .dO 

else 

dy=l.dO/dsqrt(b) 
endif 
else 

a=q ( ix+ 1 ) -q ( ix) - 1 . dO ; b=l . dO+8 . dO*a* (x-q( ix) ) 
if(b <= O.dO) then 

vkoch=l .dO 
else 

dy=l.dO/dsqrt(b) 
endif 
endif 
end function dy 
end function dp_ 
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function r(t) ! asymptote of function p(t) 
M. Cipolla, "La determinazione assintotica dell nimo numero primo 
Rend. Acad. Sci. Fis. Mat. Napoli, Ser. 3, 8 (1902), 132-166 

implicit real*8(a-h,o-z) 

common/ rcorr/ rc , drc 
r=t* (dlog(t) +dlog (dlog (t) )+ (dlog (dlog (t) ) -2 . dO) / 

dlog(t)-((dlog(dlog(t)))**2-6.d0*dlog(dlog(t))+ll.d0)/ 

(2.d0*dlog(t)**2)-l.d0)+rc 
end function r 

function dr(x) ! derivative of asymptote r(t) 
implicit real*8(a-h,o-z) 
common/ rcorr/ rc , drc 
tl = dlog(x); t2 = dlog(tl); t3 = tl**2; t4 = t3*tl 
t7 = tl-2.d0; t8 = dlog(t7) ; tl8 = t2**2; t21 = t3**2 
t35 = -56.d0-4.d0*t2*t4+4.d0*t8*tl+50.d0*tl-13.d0*t3+ 
2 . d0*t4-6 . d0*t8*t3-26 . d0*t2*t 1+4 . dO*t 18*t 1+ 
2 . d0*t2*t21+2 . d0*t8*t4-l . d0*tl8*t3+6 . d0*t2*t3+ 
28 . d0*t2-4 . d0*tl8-4 . d0*t21+2 . d0*t21*tl 
dr = 0.5d0*t35/t4/t7+drc 
end function dr 

subroutine eratosthenes (k) 
implicit real*8(a-h,o-z) 

common/ protiarithmi/qC 10000000) / rcorr/ rc , drc/ wn/wn/nn/nn 
common/ ck/vkoch 

rc=0.dO; drc=0.dO; np=0; vkoch=0.dO 
do n=2,k 
id=l 

isqrtn=dsqrt (df loat (n) ) 
id=id+l 

if (n == 2) goto 2 

if (mod (n, id) == 0) goto 3 
if (id >= isqrtn) goto 2 

goto 1 
np=np+l; q(np)=n 
continue 

enddo 

wn=df loat (np) ; nn=np; rc=p(wn)-r(wn) ; drc=dp (wn) -dr (wn) 
end subroutine eratosthenes 

subroutine primes (n) 
implicit real*8(a-h,o-z) 

common/protiarit]imi/q( 10000000) /rcorr/rc,drc/wn/wn/nn/nn 

common/ ck/ vkoch 

rc=0.dO; drc=0.dO; nn=n; wn=df loat (n) ; vkoch=0.dO 
open (222, file=' primes') 
do i=l, n/6+1 

il=6*i-5 ; i2=6*i-4 ; i3=6*i-3 ; i4=6*i-2 ; i5=6*i-l ; i6=6*i ; 
read(222, *) q(il) ,q(i2) ,q(i3) ,q(i4) ,q(i5) ,q(i6) 
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enddo 

nn=n; rc=p(wn)-r(wn) ; drc=dp (wn) -dr (wn) 
end subroutine primes 

! 

function dlogintegral(x) 

implicit real*8(a-h,o-z) 

parameter (itol=100, small=l .d-10) 

dli=0.dO; i=0 
1 i=i+l; dliw=dli; dm=(dlog(x))**i 

if (dm > small) then 

dli=dli+df loat (factorial (i-1) ) /dm 

if (dabs (dliw-dli) > small. and. i <= itol) goto 1 

endif 

dlogintegral=x*dli 
contains 

integer recursive function factorial (1) result (If) 

integer (4) 1 

lf=l 

if(l > 0) then 

If =l*f actorial(l-l) ; return 
endif 

end function factorial 
end function dlogintegral 



program test_pp_f 
implicit real*8(a-h,o-z) 

open(21, file='p.txt') ; open(22, file= 'p_.txt') 
open(23, file='dp.txt') ; open(24, file=' dp_.txt') 
n=6000 

call eratosthenes(n) ! call primes (n) 
s=5970.d0; sm=0.05d0; ns=1000 

!uses function primes(n) at n=6000; returns tab for fig 4 

s=428.d0; sm=0.01dO; ns=600 ! returns tabs for figs 1, 2 and 3 
s=154.78d0; sm=0.01dO; ns=1350 ! returns tabs for figs 6,7,8 and 9 
s=l.dO; sm=0.1dO; ns=250; ! returns tab for fig 5 

ss=s 

do i=l, ns 

s=s+sm 
write(21,*) s, p(s) 
write(23,*) s, dp(s) 
enddo 

sm= (p (s) -p (ss) ) /df loat (ns) ; s=p (ss) 

do i=l, ns 

s=s+sm 

write (22,*) s, p_(s) ! p_newt(s) 

write (24,*) s, dp_(s) 
enddo 
end 
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Application 2. 

Solution of the equation + x(2)^ = x(3)^ + 1 over primes 



! user module to the main program af xy 

subroutine FXY(m,n,np,neq,f ,x,pp,df ,yr) 
implicit real*8(a-h, o-z) 
DIMENSION X(l) ,pp(l) ,DF(1) ,YR(1) 

COMMDN/LSMH/LSMH , MQH , NSOLH/BXH/Dl , D2 , BL (600) , BR(600) 
COMMON/RETFH/LFl , LF2 , LF3 , NDAT/FR/FR/SLMH/SSVH , S3H 
go to (1,2) , np 

I run-time — section 

2 pi=dacos(-l .DO) 
goto (21, 22) , neq 

! .first. .equation 

21 f=x(l)**2+x(2)**2-x(3)**2-l.d0 ! diophantine equation 
df(l)= 2.d0*x(l); df(2)= 2.d0*x(2); df (3)=-2.d0*x(3) ! derivatives 
return 

! .second. .equation 

22 f=0.d0 
do i=l,3 

df (i)=pi*dsin(2 . dO*pi*p_(x(i) ) ) *dp_(x(i) ) ! derivatives 

f=f+(dsin(pi*p_(x(i))))**2 ! f rom-reals-to-primes equation 

enddo 

return 

! pre-execution — section — set — the — controls — to — af xy-program 

1 nn=10000; call eratosthenes (nn) ! prime number table creation 
n=3 ! number of unknowns 

m=2 ! number of equations 

np=-3 ! produce autoregularized Gauss-Newton process 
f=l.D-34 ! accuracy level for the residual f(x)-y 
lsmh=l ! Gene Golub's SVD-method 

ssvh=l.d-16 ! minimal characterist numb in SVD-method 
mqh=l ! Jorge More's addaptive scaling 

nsolh=20 ! limit of the sought solutions 

do i=l, n 

x(i)=20.d0 ! guesses 

bl(i)=2.d0; br(i)=100.d0 ! constraints (definition domain Df) 
enddo 
If 1=1; If 2=2 
return 

end 

! 

FOUND SOLUTIONS: 
Solution # 1 (Plan 2; xO # 2; eps0=2 . 88D-08) : 

K= 32 

(Df)'(fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

2.66223D-21 6.39687D-15 4.09199D-29 1.3172D+04 3.5976D-23 O.OOOD+00 4.8D-06 
Unknowns : 

X( 1)= 2.8999999979D+01 X( 2)= 2 . 3000000010D+01 X( 3)= 3. 6999999989D+01 
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Solution # 2 (Plan 2; xO # 2; eps0=2 . 88D-03) : 

K= 26 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

5.32001D-18 9.98342D-13 9.96687D-25 1.6260D+04 1.2170D-13 O.OOOD+00 4.8D-02 

Unknowns : 

X( 1)= 2.9000000270D+01 X( 2)= 2.8999999962D+01 X( 3)= 4. 1000000164D+01 
Solution # 3 (Plan 2; xO # 2; eps0=2 . 88D-01) : 

K= 50 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

7.48702D-21 1.30119D-14 1.69311D-28 2.6924D+04 7.0632D-23 O.OOOD+00 5.3D-01 

Unknowns : 

X( 1)= 4.2999999971D+01 X( 2)= 3 . 1000000004D+01 X( 3)= 5 . 2999999979D+01 
Solution # 4 (Plan 2; xO # 3; epsO=2.42D-08) : 

K= 35 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

1.74942D-20 2.45302D-14 6.01732D-28 1.3197D+04 1.3640D-13 O.OOOD+00 3.1D-05 
Unknowns : 

X( 1)= 2.2999999992D+01 X( 2)= 2.8999999964D+01 X( 3)= 3 . 6999999967D+01 
Solution # 5 (Plan 2; xO # 3; epsO=2.42D-07) : 

K= 52 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

9.75966D-22 3.70459D-15 1.37240D-29 5.0600D+03 1.7760D-14 O.OOOD+00 3.1D-04 
Unknowns : 

X( 1)= 1.9000000009D+01 X( 2)= 1 .3000000011D+01 X( 3)= 2 . 3000000013D+01 
Solution # 6 (Plan 2; xO # 3; epsO=2.42D-04) : 

K= 32 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

6.32452D-21 1 . 18087D-14 1.39445D-28 2.1246D+04 6.7279D-23 O.OOOD+00 3.1D+00 

Unknowns : 

X( 1)= 2.9000000000D+01 X( 2)= 3.6999999973D+01 X( 3)= 4.6999999979D+01 
Solution # 7 (Plan 2; xO # 4; eps0=2 . OlD-06) : 

K= 52 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

9.97306D-21 1.59910D-14 2.55714D-28 2.0899D+04 1.0602D-22 O.OOOD+00 3.1D-05 

Unknowns : 

X( 1)= 2.3000000008D+01 X( 2)= 4 . 0999999968D+01 X( 3)= 4.6999999976D+01 
Solution # 8 (Plan 2; xO # 4; eps0=2 . OlD-05) : 

K= 56 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

2.27957D-20 3.09817D-14 9.59867D-28 1.6685D+04 2.5783D-13 O.OOOD+00 3.1D-03 
Unknowns : 

X( 1)= 1.3000000029D+01 X( 2)= 4. 1000000030D+01 X( 3)= 4.3000000037D+01 
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Solution # 9 (Plan 2; xO # 5; epsO=l . 74D-08) : 

K= 40 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

2.16010D-16 1.23358D-11 1.52172D-22 2.6935D+04 1.5541D-13 O.OOOD+00 1 . lD-06 
Unknowns : 

X( 1)= 3.0999999926D+01 X( 2)= 4.3000000887D+01 X( 3)= 5 . 3000000676D+01 
Solution # 10 (Plan 2; xO # 5; epsO=l . 74D-04) : 

K= 37 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

4.52579D-21 1.02193D-14 1.04433D-28 8.8041D+03 2.8456D-14 O.OOOD+00 1 . lD-02 
Unknowns : 

X( 1)= 2.9000000019D+01 X( 2)= 1 . 1000000013D+01 X( 3)= 3 . 1000000022D+01 
Solution # 11 (Plan 2; xO # 5; epsO=l . 74D+00) : 

K= 29 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

4.96913D-21 1.02836D-14 1.05752D-28 2.0492D+04 2.3436D-13 O.OOOD+00 2.2D+01 
Unknowns : 

X( 1)= 4.3000000024D+01 X( 2)= 1 .8999999989D+01 X( 3)= 4.7000000018D+01 
Solution # 12 (Plan 2; xO # 6; epsO=2.43D-03) : 

K= 44 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

1.26391D-23 1.87462D-16 3.51422D-32 5.0621D+03 3.5534D-14 O.OOOD+00 1.3D+01 
Unknowns : 

X( 1)= 1.3000000000D+01 X( 2)= 1 . 9000000003D+01 X( 3)= 2 . 3000000003D+01 
Solution # 13 (Plan 2; xO # 7; eps0=2 . 70D-08) : 

K= 58 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

7.20035D-21 1.26460D-14 1.59923D-28 4.3148D+04 4.7940D-16 O.OOOD+00 6.7D-03 

Unknowns : 

X( 1)= 4.1000000029D+01 X( 2)= 5.2999999979D+01 X( 3)= 6 . 7000000001D+01 
Solution # 14 (Plan 2; xO # 7; eps0=2 . 70D-02) : 

K= 44 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

7.26781D-17 6.74038D-12 4.54328D-23 6.5404D+04 2.1514D-13 O.OOOD+00 6.7D-01 

Unknowns : 

X( 1)= 7.1000000315D+01 X( 2)= 4.3000000534D+01 X( 3)= 8 . 3000000546D+01 
Solution # 15 (Plan 2; xO # 9; eps0=2 . 73D-03) : 

K= 28 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

2.31418D-21 6.77216D-15 4.58621D-29 8.8041D+03 3.1949D-14 O.OOOD+00 1.6D+00 
Unknowns : 

X( 1)= 1.1000000015D+01 X( 2)= 2 . 9000000013D+01 X( 3)= 3 . 1000000017D+01 
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Solution # 16 (Plan 2; xO # 10; eps0=2 . 78D-05) : 

K= 71 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

1.32835D-20 2.01058D-14 4.04243D-28 4.3148D+04 9.9131D-23 O.OOOD+00 1.7D-04 

Unknowns : 

X( 1)= 5.2999999967D+01 X( 2)= 4. 0999999995D+01 X( 3)= 6 . 6999999970D+01 
Solution # 17 (Plan 2; xO # 13; eps0=2 . 58D+00) : 

K= 36 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

7.87327D-20 5.71405D-14 3.26504D-27 7.4410D+04 2.8416D-13 O.OOOD+00 8.1D+01 

Unknowns : 

X( 1)= 7.9000000002D+01 X( 2)= 4 . 0999999930D+01 X( 3)= 8 . 8999999970D+01 
Solution # 18 (Plan 2; xO # 15; epsO=2.47D-03) : 

K= 55 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

1.70828D-20 2.50831D-14 6.29162D-28 4.7012D+04 5.1160D-13 O.OOOD+00 3.4D+00 
Unknowns : 

X( 1)= 7.1000000031D+01 X( 2)= 1 .7000000020D+01 X( 3)= 7 . 3000000035D+01 
Solution # 19 (Plan 2; xO # 16; eps0=2 . 29D-09) : 

K= 47 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

6.22186D-21 1.25141D-14 1.56603D-28 1.6703D+04 1.2445D-13 O.OOOD+00 1.4D-05 
Unknowns : 

X( 1)= 4.1000000024D+01 X( 2)= 1 .3000000008D+01 X( 3)= 4.3000000025D+01 
Solution # 20 (Plan 2; xO # 21; epsO=l . 97D-08) : 

K= 34 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

1.02797D-20 1.75525D-14 3.08090D-28 2.1356D+04 1.5072D-13 O.OOOD+00 2.1D-05 

Unknowns : 

X( 1)= 4.0999999971D+01 X( 2)= 2.2999999991D+01 X( 3)= 4.6999999971D+01 

Singular values of the matrix Z=(Df ) 'Df at the best approximation X( 34) 
(errSVD= 0) : 

1) O.OOOOOOOOOOOOOD+00 2) 2 . 0194839173658D-28 3) 1 .8088976985834D+04 

Cond(Z(x( 34))=104 ; GMcond(Z(x( 34)))= 45 ; nullity(Z(x( 34)))= 2 . 

Residual_tab of the vector r=y-f(x( 34)) (M= 2): 
O.OD+00, 1; -1.8D-14, 1; 

Final report : 202 iterative processes were tried, 
9093 times subroutine FXY was called, 
8911 iterations were produced, and 
20 solutions were found. 
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Application 3. 

Solution of the system x(l)^ + x(2)^ = x(3)^ + 1, x(3) — x(l) = 2 over twins 

! user ' s module — to — the — main — program — af xy 

subroutine FXY(m,n,np,neq,f ,x,pp,df ,yr) 
implicit real*8(a-h, o-z) 

DIMENSION X(l) ,pp(l) ,DF(1) ,YR(1) 

COMMON/LSMH/LSMH , MQH , NSOLH/BXH/Dl , D2 , BL (600) , BR(600) 
COMMON/RETFH/LFl , LF2 , LF3 , NDAT/SLMH/SSVH , S3H 
go to (1,2) , np 

I run-time — section 

2 continue; pi=dacos(-l .DO) 
goto (21, 22, 23) , neq 
! . . .first. .equation 

21 f=x(l)**2+x(2)**2-x(3)**2-l.d0 ! diophantine equation 
df(l)= 2.d0*x(l); df(2)= 2.d0*x(2); df (3)=-2.d0*x(3) ! derivatives 
return 

! . . .second. .equation 

22 f=x(3)-x(l)-2.d0 ! twin-pair equation 
df(l)=-l.dO; df(2)= O.dO; df(3)= l.dO ! derivatives 
return 

! . . .third. .equation 

23 f=0.dO 
do i=l,3 

df (i)=pi*dsin(2 . dO*pi*p_(x(i) ) ) *dp_(x(i) ) ! derivatives 

f=f+(dsin(pi*p_(x(i))))**2 ! f rom-reals-to-primes equation 
enddo 

return 

! pre-execution — section — set — the — controls — to — af xy-program 

1 continue 

nn=10000; call eratosthenes(nn) ! prime number table creation 
n=3 ! number of unknowns 

m=3 ! niimber of equations 

np=-3 ! produce autoregularized Gauss-Newton process 
f=l.D-34 ! accuracy level for the residual f(x)-y 

lsmh=l ! Gene Golub's SVD-method 

ssvh=l.d-16 ! minimal characterist numb in SVD-method 
mqh=l ! Jorge More's addaptive scaling 

nsolh=10 ! limit of the sought solutions 

do i=l, n 

x(i)=50.d0 ! guesses 

bl(i)=2.d0; br(i)=i01.dO ! constraints (definition domain Df) 
enddo 
If 1=1; If 2=2 
return 

end 
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FOUND SOLUTIONS: 
Solution # 1 (Plan 2; xO # 1; eps0=3 . OOD-02) : 

K= 31 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

9.94498D-20 7.96082D-14 6.33746D-27 4.7014D+04 6.3947D-14 O.OOOD+00 3.6D-01 

Unknowns : 

X( 1)= 7.0999999937D+01 X( 2)= 1 . 6999999993D+01 X( 3)= 7 . 2999999937D+01 
Solution # 2 (Plan 2; xO # 1; eps0=3 . OOD+00) : 

K= 36 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

1.84559D-21 5.58920D-15 3.12392D-29 1.6686D+04 7.1257D-15 O.OOOD+00 7.2D+00 
Unknowns : 

X( 1)= 4.1000000017D+01 X( 2)= 1 .3000000003D+01 X( 3)= 4.3000000017D+01 
Solution # 3 (Plan 2; xO # 2; eps0=2 . 88D-08) : 

K= 34 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

9.06718D-24 1.65074D-16 2.72495D-32 4.7800D+02 1.7750D-15 O.OOOD+00 3.1D-04 
Unknowns : 

X( 1)= 5.0000000028D+00 X( 2)= 5 . OOOOOOOOllD+00 X( 3)= 7 . 0000000028D+00 
Solution # 4 (Plan 2; xO # 2; eps0=2 . 88D-07) : 

K= 34 

(Df)'(fx-y) fx-y chi.sqr. (Df)'Df l.p.prec. eps en/r 

1.14803D-22 8.85711D-16 7.84485D-31 1.6145D+03 3.5507D-15 O.OOOD+00 3.1D-04 

Unknowns : 

X( 1)= 1.0999999993D+01 X( 2)= 6.9999999981D+00 X( 3)= 1 . 2999999993D+01 
Solution # 5 (Plan 2; xO # 2; eps0=2 . 88D-03) : 

K= 33 

(Df ) ' (fx-y) fx-y chi.sqr. (Df ) 'Df l.p.prec. eps en/r 

4.85228D-22 2.29735D-15 5.27781D-30 8.8060D+03 3.5523D-15 O.OOOD+00 3.1D-01 

Unknowns : 

X( 1)= 2.9000000011D+01 X( 2)= 1 . 1000000002D+01 X( 3)= 3 . lOOOOOOOllD+01 

Final report : 300 iterative processes were tried, 
6108 times subroutine FXY was called, 
6027 iterations were produced, and 
5 solutions were found. 
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